% B1, B2 & C1, C2

function [coef, p1] = func_coef_of_fields(m, CONSTS)

    debug_output = false;

    p_low = 1 + 1e-10;
    p_high = 10;
    use_combination = false;
    num_intervals = 200000; % 2000
    
    f = @(p)imag(func_delta(p, m, CONSTS));
%     f1 = @(p)real(func_delta(p, m, CONSTS));

    val_low = f(p_low);
    val_high = f(p_high);
    
%     if (val_low * val_high > 0)
%         coef = [0; 0; 0; 0];
%         p1 = NaN;
%         return;
%     end

%     z = [fzeroall(f, [p_low, p_high], num_intervals); ...
%         fzeroall(f1, [p_low, p_high], num_intervals)];

    z = fzeroall(f, [p_low, p_high], num_intervals); 


    if(debug_output)
        z_det = fzeroall(@(p)imag(det(matrix_of_boundary_conditions(p, m, CONSTS, use_combination))), ...
            [1, p_high], num_intervals, false);
        fprintf('Analytical zero = %f; Numerical zero = %f\n', z(1), z_det(1));
    end
    
    if (size(z, 1) <= 0)
        coef = [0; 0; 0; 0];
        p1 = NaN;
        return;
    end
    
    p1 = z(1); % was 1
    mat = matrix_of_boundary_conditions(p1, m, CONSTS, use_combination);
    det(mat);
    rank(mat);
    K = mat(1:3, 1:3);
    det(K);
    rank(K);
    C2 = (1/1i)*det(K);
    f = - C2 * mat(1:3, 4);
    u = K^(-1) * f;
    B1 = u(1);
    B2 = u(2);
    C1 = u(3);
    %C2
    
    transp = false;
    
    if(transp)
        
        z_T = fzeroall(@(p)imag(det(matrix_of_boundary_conditions_T(p, m, CONSTS, use_combination))), ...
                           [1, p_high], num_intervals);
        p1_T = z_T(1);
        mat_T = matrix_of_boundary_conditions_T(p1_T, m, CONSTS, use_combination);
        det(mat_T);
        rank(mat_T);
        K_T = mat_T(1:3, 1:3);
        det(K_T);
        rank(K_T);
        C2_T = (1/1i)*det(K_T);
        f_T = - C2_T * mat_T(1:3, 4);
        u_T = K_T^(-1) * f_T;
        
        fprintf('B1 = %e + %ei; B1_T = %e + %ei\n', real(B1), imag(B1), real(u_T(1)), imag(u_T(1)));
        fprintf('B2 = %e + %ei; B2_T = %e + %ei\n', real(B2), imag(B2), real(u_T(2)), imag(u_T(2)));
        fprintf('C1 = %e + %ei; C1_T = %e + %ei\n', real(C1), imag(C1), real(u_T(3)), imag(u_T(3)));
        fprintf('C2 = %e + %ei; C2_T = %e + %ei\n', real(C2), imag(C2), real(C2_T), imag(C2_T));
        
    end

    epsilon = 1e-10;
    %coef = [B1; B2; real(C1); real(C2)];
    %coef = clear_small_values([B2; B1; C1; C2], epsilon);
    coef = [B1; B2; C1; C2];

    eq = clear_small_values(mat * coef, epsilon);
    
    if(debug_output)
        fprintf('equations for p1  %e \n', eq);
    end

    if(debug_output)
        
        E_z_inside_a  = mat(2, 1)*B1 + mat(2, 2)*B2;
        E_z_outside_a = mat(2, 3)*C1 + mat(2, 4)*C2;

        H_z_inside_a  = mat(1, 1)*B1 + mat(1, 2)*B2;
        H_z_outside_a = mat(1, 3)*C1 + mat(1, 4)*C2;

        E_phi_inside_a  = mat(3, 1)*B1 + mat(3, 2)*B2;
        E_phi_outside_a = mat(3, 3)*C1 + mat(3, 4)*C2;

        H_phi_inside_a  = mat(4, 1)*B1 + mat(4, 2)*B2;
        H_phi_outside_a = mat(4, 3)*C1 + mat(4, 4)*C2;

        fprintf('r=a: E_z_inside   = %e + %ei      E_z_outside   = %e + %ei \n',  ...
                               real(E_z_inside_a),  imag(E_z_inside_a), ...
                               real(E_z_outside_a), imag(E_z_outside_a));
        fprintf('r=a: H_z_inside   = %e + %ei       H_z_outside   = %e + %ei \n',  ...
                               real(H_z_inside_a),  imag(H_z_inside_a), ...
                               real(H_z_outside_a), imag(H_z_outside_a));
        fprintf('r=a: E_phi_inside = %e + %ei       E_phi_outside = %e + %ei \n',  ...
                               real(E_phi_inside_a),  imag(E_phi_inside_a), ...
                               real(E_phi_outside_a), imag(E_phi_outside_a));
        fprintf('r=a: H_phi_inside = %e + %ei      H_phi_outside = %e + %ei \n',  ...
                               real(H_phi_inside_a),  imag(H_phi_inside_a), ...
                               real(H_phi_outside_a), imag(H_phi_outside_a));
    end
    
end